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We describe an implementation of Hedin's GW approximation for molecules and clusters, the 
complexity of which scales as 0(N 3 ) with the number of atoms. Our method is guided by two 
I strategies: i) to respect the locality of the underlying electronic interactions and ii) to avoid the 

singularities of Green's functions by manipulating, instead, their spectral functions using fast Fourier 
(**-) transform methods. To take into account the locality of the electronic interactions, we use a local 

basis of atomic orbitals and, also, a local basis in the space of their products. We further com- 
press the screened Coulomb interaction into a space of lower dimensions for speed and to reduce 
memory requirements. The improved scaling of our method with respect to most of the published 
methodologies should facilitate GW calculations for large systems. Our implementation is intended 
as a step forward towards the goal of predicting, prior to their synthesis, the ionization energies and 
electron affinities of the large molecules that serve as constituents of organic semiconductors. 
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I. INTRODUCTION 



Lars Hedin's GW method [T] is an approximate treatment of the propagation of electrons in condensed matter where 
an electron interacts with itself via a Coulomb interaction that is screened by virtual electron-hole pairs. In periodic 
semiconductors, the GW approximation is known to lead to surprisingly accurate gaps [2], while for finite clusters 
and molecules it provides qualitatively correct values of ionization energies and electron affinities [3J. Hedin's GW 
approximation is also needed, as a first step, when using the Bethe-Salpeter equation to find the optical properties of 
systems in which the Coulomb interaction is only weakly screened. 

The present work is motivated by the rapid progress, during the last decade, in the field of organic semiconductors, 
especially in organic photovoltaics and organic luminescent diodes [3]. To optimize such systems, it would be useful 
to know the key properties of their molecular constituents before actually synthesizing them. In order to make such 
predictions it is necessary to develop algorithms with a favorable complexity scaling, since many of the technologically 
relevant molecules are fairly large. The method presented here is a step forward along this direction. Its 0(N 3 ) 
scaling, with N the number of atoms in the molecule, is an improvement over most existing methodologies. 

While computational techniques for treating the GW approximation for clusters and molecules have become sophis- 
ticated enough for treating molecules of interest in photovoltaics [5] or in the physiology of vision [5] , such calculations 
remain computationally expensive. The scaling with the number of atoms of these recent calculations has not been 
published. However, in many cases it is unlikely to be better than 0(N 4 ) [7]. A recently published method for com- 
puting total energies of molecules that uses the random phase approximation also has 0(N 4 ) scaling [5j. Actually, at 
this point it is difficult for us to envisage a scaling exponent less than three because the construction of the screened 
Coulomb interaction — the central element of the GW approach — requires inverting a matrix of size 0(N) which, 
in general, takes 0(N 3 ) operations. 

The algorithm described in this paper is based on two main ingredients: i) respecting the locality of the underlying 
interactions and, ii) the use of spectral functions to describe the frequency/time dependence of the correlators. The 
latter ingredient allows for the use of the fast Fourier transform (FFT) to accelerate the calculations, while the former 
idea of respecting locality has been also at the heart of other efficient GW methods, like the successful "space-time 
approach" for periodic systems {[)■. 

Our method is based upon the use of spatially localized basis sets to describe the electronic states within the 
linear combination of atomic orbitals (LCAO) technique. In particular, we have implemented our method as a post 
processing tool of the SIESTA code |10) . although interfaces with other LCAO codes should be simple to construct. 
The precision of the LCAO approach is difficult to control and to improve, but a basis of atom centered local orbitals is 
useful for systems that are too large to be treated by plane- wave methods [TT]. In order to solve Hedin's equations we 
construct a basis set that gets rid of the over completeness of the orbital products while keeping locality. In molecular 
computations this is frequently done through a fitting procedure (using Gaussians or other localized functions). We 
use an alternative mathematical procedure [121 113) that dispenses with this fitting and defines a basis of dominant 
products. The basis of dominant products was instrumental to develop an efficient linear response code for molecular 
absorption spectra |14[ 115] . In the present paper we have developed an additional, non local compression technique 
that further reduces the size of the product basis. The compression allows to store the whole matrix representation of 
the screened Coulomb interaction at all times/frequencies and needs much less memory. Moreover, the compression 
strongly accelerates the calculation of the screened Coulomb interaction because it involves a matrix inversion. This 
leads to a gain in computational efficiency which is even more important than that associated with the reduction of 
the needed memory. 

Of course there are other methods that use a localized basis different from LCAO and, thus, equally appropriate for 
dealing with clusters and molecules while exploiting locality. One method uses a lattice in real space [16 . Another 
method uses wavelets that represent a useful compromise between localized and extended states [17]. Localized 
Wannier orbitals obtained from transforming plane waves [TH] have also been used in GW calculations [TH1 [20] • In 
this paper we use a basis of dominant functions to span the space of products of atomic orbitals [12] and we use a 
compression scheme to deal with the screened interaction. It is clear, however, that some of the ideas and techniques 
of the present paper can be combined with the alternative approaches quoted above. 

The actual implementation of the algorithm that we report in this paper can be considered as a "proof of principle" 
only and the prefactor of our implementation leaves room for further improvement. Therefore, we validate our method 
with molecules of moderate size (we consider molecules of only up to three aromatic rings: benzene, naphthalene and 
anthracene) , leaving further improvements and applications to molecules of larger sizes for a future publication. 

This paper is organized as follows. In section ITT] we recall the equations of the GW approximation. In section |III| 



we rewrite the GW approximation for molecules in tensorial form. Section IV describes the instantaneous component 
of the self-energy, while in section [V] we describe a spectral function technique for solving these tensorial equations. 
Section [VT| describes our GW results for benzene, a typical small molecule. Section |VTTJ describes our algorithm for the 
compression of the screened Coulomb interaction that is needed to treat larger molecules. Section VIII explains how 
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0(N 3 ) scaling can be maintained for large molecules by alternatively compressing and decompressing the Coulomb 
interaction. Section IX presents a summary of the entire algorithm for performing GW calculations. In section [X] we 
test our method on naphthalene and anthracene, and our conclusions are presented in section |XI) 



II. ELEMENTARY ASPECTS OF HEDIN'S GW APPROXIMATION 



The one-electron Green's function of a many-body system has proved to be a very useful concept in condensed 
matter theory. It allows to compute the total energy, the electronic density and other quantities arising from one- 
particle operators. The one-electron Green's function G(r, r', t) has twice as many spatial arguments as the electronic 
density but it remains a far less complex object than the many-body wave function. Furthermore, Hedin has found an 
exact set of equations for a finite set of correlation functions of which the one-electron Green's function is the simplest 
element. This set of equations has not been solved so far for any system whatsoever. However, as a zeroth order 
starting point to his coupled equations, Hedin suggested the very successful GW approximation for the self-energy 
T,(r,r',t). This approximation describes the change of the non interacting electron propagator Go(r,r',t) due to 
interactions among the electrons. With the help of a self-energy, one can find the interacting Green's function from 
Dyson's equation 

G = Go + GoEGo + GoSGoSGo + . . . = — — [ — — , (1) 

G - L 

where products and inversions in the equation must be understood in an operator sense as required in many-body 
perturbation theory. 

In Hcdin's GW approximation, the interaction of electrons with themselves is taken into account by the following 
self-energy 



Z(r,r',t)=iG (r,r',t)W(r,r',t), (2) 

where W(r,r',t) is a screened Coulomb interaction. 

The key idea of Hedin's GW approximation [1] is to incorporate the screening of the Coulomb interaction from the 
very beginning in a zeroth order approximation. Let v(r,r') = \r — r'| _1 be the bare Coulomb interaction and let 

Xo( r , r ' , t — f) = svlr 1 ' 1 }') k e the density response 6n(r, t) of non interacting electrons with respect to a change of the 
external potential 6V(r' ,t'). Hedin then replaces the original Coulomb interaction v(r,r') by the screened Coulomb 
interaction W(r,r',u>) within the random phase approximation (RPA) |2l] 

W(r, r', u) = T7 : ~ f — — ■ — tt — ; Mr'", r'), (3) 

v ' ' ' d(r -r"')-v(r,r")xa(r" ' ,r"> \u) v ; yj 

here and in the following we assume integration over repeated spatial coordinates (in our case r" and r ) on the 
right hand side of an equation if they do not appear on its left hand side. This convention makes our equations more 
transparent without introducing ambiguities and it is analogous to the familiar Einstein's convention of summing over 
repeated indices. 

We can justify the expression (|3J) by considering an internal screening field SVi n duced( r ■> w) that is generated by an 
extra external field SV cx t crna i(r , oj) 

SV tota i(r,0j) = ^external (r,Ul) + <5Vj nd uccd (r, Cj) , (4) 

where 

SV induccd (r,uj) = v(r 7 r")Sn induccd (r" ,u) = v(r, r")xa(r", r'" , oS)5V to%a \{r"' , w). 
As a consequence we obtain a frequency dependent change of the total potential 

5VW(T,w) = Tf -, 1 /n -, „ Jf. r ^external (t-"',U>). (5) 

o(i — r'") — v{r, r")xo{ r > r , 

If we assume that large fields are screened the same way as small field changes, then we may replace 5V ex teinai(f , uj) 
by the singular Coulomb interaction v(r,r') an d we obtain the screened counterpart W[r,r',uS) of the original bare 
Coulomb interaction as in eq (|3|. 
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Because of the relation 

i X o(r, r', t) = 2G (r, r', t)G (r', r, -t) (6) 

the screening in eq (|3| may be interpreted as being due to the creation of virtual electron-hole pairs. The screening 
by virtual electron-hole pairs is the quantum analogue of classical Debye screening in polarizable media |22j . The 
factor of 2 in eq ([6| takes into account the summation over spins. 

Many body theory uses Feynman-Dyson perturbation theory |23j and the latter is formulated in terms of time 
ordered correlators. For instance, a Green's function is represented as a time ordered correlator of electron creation 
ip + (r,t) and annihilation ip(r,t) operators 

iG(r,r',t-i') = ^(^-i')(0|V'(r,t)^ + (r',t')|0)-^'-t)(0|^ + (r',t')VXr-^)|0), (7) 

where the minus sign is due to Fermi statistics, |0) denotes the electronic ground state, and where 6(t) denotes 
Heaviside's step function. This completes our formal description of Hedin's GW approximation. 

In practice, Hedin's equations are solved "on top" of a density functional or Hartree-Fock calculation. The framework 
of density functional theory (DFT) [2U [25] already includes electron correlations at the mean-field level via the 
exchange correlation energy E xc [n(r)], where [n] denotes the functional dependence of E xc on the electron density. 
DFT calculations are usually performed using the Kohn-Sham scheme [25] . in which electrons move as independent 
particles in an effective potential. The Kohn-Sham Hamiltonian i?KS reads 

H KS - -\^ 2 + V KS , (8) 

8E 

Vks = Vest + VHartrco + V xc , where V xc (r) 



Sn(r) 



To avoid including the interaction twice, the exchange correlation potential V xc (r) must be subtracted from E(r, r' , t) 
in eq ^ when using the output of a DFT calculation as input for a GW calculation. This is done by making the 
replacement 

E(r, r', t) -> E(r, r', t) - S(r - r')S(t)V KC (r) (9) 

in Dyson's equation ([I]). 

Our aim is to compute the electronic density of states (DOS) that is defined as the trace of the imaginary part of 
the electron propagator 

p(oj + ie) = -^Im J G(lj + ie, r, r)d 3 r. (10) 

The electronic DOS p(u> + ie) can be compared with experimental data from direct and inverse photo-emission 26 . 
From it, we can read off the energy position of the highest occupied and the lowest unoccupied molecular orbitals 
(HOMO and LUMO) or, alternatively, the ionization energy and the electron affinity. 
Finally, let us list the equations that define the GW approximation 



iXo( T ', r ' 1 1) = 2Go(r, r', t)Go(r', r, — t); free electron response 

W{r, r', oj) = [S(r - r'") - v(r, r")x [r", r'", a;)]" 1 v(r"', r'); RPA screening 

E(r,r',t) = iG (r,r',t)W(r,r',t); GW self-energy 

G _1 (r, r' , u) + ie) = Gq 1 (r, r', oj + ie) — S(r, r' , oj + ie). Dyson equation 



(11) 



The next two sections describe the tensor form of equations (11) as well as the main ingredients of our imple- 
mentation of the GW approximation as embodied in equations (11) for the case of small molecules. Later sections 
will describe the compression/decompression of the Coulomb interaction that is needed for treating large molecules 
without over flooding the computer memory. 



III. TENSOR FORM OF HEDIN'S EQUATIONS 



In order to compute the non interacting Green's function (7]), we will use the LCAO method where one expresses 
the electron operator in terms of a set of fermions c a (t) that belong to localized atomic orbitals [57] 



^(r,t)^^2r(r)c a (t). 



(12) 
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Such a parametrization is parsimonious in the number of degrees of freedom, although its quality is difficult to control 
and to improve in a systematic way. 

The output of a DFT calculation (that serves as input for the GW calculation) is the Kohn-Sham Hamiltonian 
H ab and the overlap matrix S of the LCAO basis functions f a (r) [25]. One may use the eigenvectors {X^} of the 
Kohn-Sham Hamiltonian 

H ab X b E = ES ab X^ (13) 
to express the (time-ordered) propagation of electrons between localized atomic orbitals 



yfi vE 

G° fc (-±k)=E " 



w ± is - E 



(14) 



In this paper we measure energies relative to a Fermi energy, so that E < (E > 0) refers to occupied (empty) 
states, respectively, and the infinitesimal constant e shifts the poles of the Greens function away from the real axis. 
Moreover, to avoid cluttering up the notation, we will often use Einstein's convention of summing over repeated 
indices, as in eq (13). 



The set of equations ( 11 1 contains correlation functions, such as the density response function x(r, r' , t) that must 



be represented in a basis of products of atomic orbitals 



(15) 



\x{r,r',t-t') = 9{t - t')(0\n(r,t)n(r',t')\0) + 9(t' - t)(0\n(r' ,t')n(r,t)\0). 
Indeed, by virtue of eq (12 1, the electronic density n(r,t) — ip + (r,t)ijj(r,t) involves products of atomic orbitals 

n{r,t) = ff> + (r,t)tl>(r,t) = £ f a (r)f b (r)c+(t)c b (t). 

a, 6 

The set of products {f a (r)f b (r)} is well known to be strongly linearly dependent [25]. As an improved solution 
of this very old technical difficulty, we previously developed an algorithm to construct a local basis of "dominant 
products" F^(r) that i) spans the space of orbital products with exponential accuracy and which it) respects the 
locality of the original atomic orbitals [12]. Moreover, the products of atomic orbitals f a (r)f b (r) relate to dominant 
products F^(r) via a product vertex V? 



f a (r)f b (r) = J2v, ab F^r). 



(16) 



Because the dominant products F^(r) are themselves special linear combinations of the original products, arbi- 
trary extra fitting functions do not enter into this scheme. In order to respect the principle of locality, the above 
decomposition is carried out separately for each pair of atoms, the orbitals of which overlap. By their construction, 
the set of coefficients V? is sparse in the sense that V^ b ^ only if a,b,/i all reside on the same atom pair |12j . 
In the construction of the dominant product basis, we made use of Talman's algorithms and computer codes for the 
expansion of products of orbitals about an arbitrary center and we also used his fast Bessel transform 30J. 

To rewrite the defining equations of the GW approximation (111 in our basis, we expand both G(r,r',t — t') and 
£(r,r',i - if) in atomic orbitals f a (r) [13] 



G{ry i t-t') = G ab (t-t')t{r)f b {r 1 )- 
X(ry,t-t') = Z ab (t-t')r(r)f b (r'). 
We also develop the screened Coulomb interaction in dominant products 

W^(t-t')= [ d z rd z r'F^{r)W{r,r',t^t')F v {r l ). 



(17) 



(18) 



Using eqs (16 17 18) it is easy to show [14] that Hedin's equations take the following tensorial form in our basis 

raa'r<a u^rbb'nQ , + v free electron response (19) 

RPA screening (20) 



i X %,(t) = 2V™'G° ab (t)V bb 'G Q a , b ,(-t); 
W^(u) ' 



X ab (t)=iV« a 'G° a , b ,(t)V b ' b W^(t); 



GW approximation 
Dyson's equation 



(21) 
(22) 



FIG. 1. Feynman diagram for the GW self-energy expressed in our local LCAO and dominant products basis. 



Here v^ v denotes the Coulomb interaction v^ u = J cPrcPr' F^{r) \ r _ r ,\ F v {r') which, due to its positivity and sym- 
metry, we also refer to as a "Coulomb metric" . Indices are raised or lowered using either the overlaps of the dominant 
functions F^(r) or the overlaps of the atomic orbitals f a {r) and which are defined as follows 



d 3 rF fl (r)F l, (r), S" 



ab 



d 3 rf a (r)f b (r). 



(23) 



Figure 1 shows the Feynman diagram corresponding to eq (|2l|). The local character of the product vertex V° 
emphasized in this figure. 



is 



IV. THE INSTANTANEOUS PART OF THE SELF-ENERGY 



When the screened Coulomb interaction W^ v in eq (20) is expanded as a function of v\°, its first term is the bare 
Coulomb interaction v^ v S(t — t') and the corresponding self-energy in eq (21 ) is frequency independent. In textbook 
treatments of the theory of the electron gas, it is explained [33] that the Green's function G a b(t — f) of the electron 
gas must be defined, at t — t' , by setting t — t' = 0~ or by first annihilating and then creating electrons. Using this 



to the exchange operator 



prescription and eqs (20 21 ) one finds the following result for the frequency-independent self-energy that corresponds 



In the frequency domain, the last operator becomes a frequency independent matrix 



Ef = V™ 



£ x5x§v*** 

E<0 



(24) 



which can be computed in 0(N 3 ) operations by using the sparsity of the product vertex V® a . 

For small molecules and clusters, the instantaneous self-energy that incorporates the effects of electron exchange 
may dominate over the remaining frequency dependent self-energy. If this is the case, we may substitute E° & into eq 
(21 ) and finish the calculation by computing the DOS 



p(u + ie) = — S ab lmG ba {cj + is), 

where we have emphasized the non-orthogonality of the basis orbitals by the explicit inclusion of the overlap S ab . 
However, the frequency dependent part of the self-energy contains correlation effects that significantly improve the 



calculation quantitatively and qualitatively as we demonstrate in sections VI and |XJ Therefore, we shall present our 
approach for the frequency dependent part of the self-energy in the next section. 



USING SPECTRAL FUNCTIONS TO COMPUTE THE SELF-ENERGY 



One might want to solve the eqs p0] ) and ( 22 ) directly as matrix valued equations in time t and to use FFTs [3T] 
to shuttle back and forth between the time and frequency domains. Unfortunately, however, this direct approach is 
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doomed to fail — the functions {G a b(t), T, ab (t), W^ v {t)} are too singular at t — to be multiplied together. We will 
now show how spectral functions come to the rescue and allow us to i) respect locality in our calculations and to ii) 
accelerate our calculation by means of FFT. 

Let us consider the energy dependent density matrix 

Pab (u) = J2^X^S(u J ~E) (25) 



and rewrite the electronic propagator ( 14 ) with the help of it 



Gl b {u± l e)= ^ ds P " b(s) . (26) 

Integral representations such as these are very useful, even in finite systems where the spectral weight is concentrated 
at isolated frequencies [131 [H]. Because a spectral function is broadened by the experimental resolution e, it can 
be represented on a discrete mesh of frequencies, with the distance between mesh points somewhat smaller than e. 
All the response functions considered in the present paper have a spectral representation because their retarded and 
advanced parts taken together define a single analytic function in the complex frequency plane with a cut on the real 
axis. 

A spectral representation is merely a rather thinly disguised Cauchy integral as we can see by considering the 
Cauchy integral representation of the electronic Greens function 

W^/^, (27) 

where C is a path surrounding the point z with Imz > in an anti-clockwise direction. If the point z — oo is regular, 
the complex plane may be treated like a sphere and we may deform the contour on this sphere in such a way that 
it wraps around the cut on the real axis in a clockwise direction. Finally, because Green's functions take mutually 
hermitian conjugate values G a b{z*) = G* ba (z) across the cut on the real axis, the above integral can the rewritten as 

G ab {z) = [ ds^^, Pab (z) = -hmG° ab (z) (28) 

J Z — S 7T 

with z on the upper branch of the cut. In writing the preceding equation we have used the simplifying feature that 
the electronic Green's function is a symmetric matrix in our real representation of angular momenta (the same is 
true for the screened Coulomb interaction) . In the following, we will always reconstruct correlation functions such as 
{G ab (uj), S ab (w), W^ v (uj)} from their imaginary part or from their spectral functions. The time ordered version of 
such correlators is determined above (below) the real axis for positive (negative) frequencies, respectively. 

A. The spectral function of a product of two correlators 

The well known convolution theorem [31 tells us that the spectral content of a product of two signals is the 
convolution of the spectral contents of its factors. The situation is entirely analogous for Green's functions and the 
other correlators considered here and their spectral functions. To see this, we use the spectral representations of the 



time ordered factors G ab {t), T, ab (t), W^ u {t) (the quantities in eqs ( 19 22 1 are time ordered) 



oo fO 

+ 



G ab (t) = -W(t) / dap+ (s)e- ls 4 +i0(-t) ds p~ b ( a ) e - ,flt ; 





S af, (i) = -i0(t) / da of {s)e- ist +i6(-t) / dsa ab {s)e- ist ; (29) 

oo 



W lv {t) = -i6{t) d S1 ^{s)e' lst + i9{-t) ds^(s)e- is \ 

JO J-oc 

where "positive" and "negative" spectral functions define the whole spectral function by means of Heaviside functions. 
For instance, the spectral function of the electronic Green's function reads 

p ab ( s ) = 6(s)p+ b (s)+9(-s)p- b (s). 
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These representations can be checked by transformin g (fo r example) the representation of G a b(t) into the frequency 
domain and by comparing with the known expression (26). We then compute S ab (t) from eq (21) and compare the 
result with the second of eqs ( 29 ) The spectral function of the self-energy is seen to have the expected convolution 
form 



a a _ b {s) 



oo poo 



JO 




5(si + s 2 - s) V^p+ hl {sx)V^ v {s 2 )d Sl ds 2 , 



6( Sl + s 2 - s)V™ /g 6 ,(si)V? b i^(s2)d Sl ds 2 



(30) 



— oo J — oo 



Note that, as commented above, the V^ a matrices are sparse and respect spatial locality. Finally, we can easily 
construct the full self-energy from its spectral functions o± b (s) by a Cauchy type integral 



oo ab 



a ab (s)ds 
uj ± is — s 



(31) 



By entirely analogous arguments, we can find the spectral function of the non interacting response Xuv from eq 
pi 



lo Jo 



OO rOO 



V^pU^W^Pdai-^)^ +s 2 - s)d Sl ds 2 , for s > 0; 



xL(w + ie) 



ds a ^ {3) , 
„ us + le — s 



for all s; 
for w > 0. 



(32) 



We implemented the convolutions in eqs (30 32) conveniently by FFT without encountering any singularities. Please 
observe that analytic continuations are not needed in our approach. 

We have seen in this subsection that the locality of the expressions for T, ab (t) and m ec t s 

rtl9h and rt2lh can be 

taken into account without multiplying singular Green's functions and by focusing instead on the spectrarfunctions 
of their products. 



The second window technique 



Although we only need results in a suitable low energy window —A < uj < X of a few electron volts, eqs (31 32) 
show that high energy processes at \uj\ > |A| influence quantities at low energies, such as, for example, the self-energy. 
Therefore, these high energy processes cannot be ignored and we need the imaginary part of the screened Coulomb 
interaction W^ v not only for small |cj| < A but also for larger frequencies. To find the imaginary part of W^ v , we 
also need, in view of eq (|2ll), the non interacting response both at small and at large frequencies. 

Let us see in the case oF the density response, how the necessary spectral information can be obtained from two 
separate calculations in two distinct frequency windows |13j . In the large spectral window —A < us < A a low 
resolution calculation with a large broadening (and, therefore, a coarse grid of frequencies) is sufficient to find at 
large energies \u\ > |A| 

X^(u + m aYgc ) = ds — — ^ . 



(33) 

l^largc S 

To get correct results in the low energy window — A < ui < A we must take into account the spectral weight in this 
window 

X^(w + l£small) = / ds — TT^ . , , , 

J-\ U + Ismail - S yj-A J\ J UJ + l£l arg0 - S (34) 

small window/ ■ • \ , [ largo window/ , • \] 

~ \ W It-smallj f [Xftu -T ItlargcJJ truncatcd spcctra i f unct ion ' 

Instead of doing the second Cauchy integral directly, we construct xjf^ sc wmdow from the large spectral window, using 
spectral data that are truncated for \s\ < A to avoid counting the spectral weight in —A < ui < A twice. Moreover, the 
broadening constant e is set differently in the first and second windows and corresponds to the spectral resolution in 
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these windows. The broadening constant e is chosen automatically in each frequency window, by setting e — 1.5Aw, 
where Aw is the distance bet ween two points on the corresponding frequency grid. We use the second window 



technique (presented in eq ([34]) for the case of the density response) again in the calculation of the self-energy S ab (w), 
where we also need the screened interaction in two windows. We combine the spectral functions of the self-energy 
in exactly the same way as for the density response. For the cases considered here, computations using two spectral 
windows were up to one order of magnitude faster than computations using a single spectral window. 

VI. TESTING OUR IMPLEMENTATION OF GW ON A SMALL MOLECULE 

The methods presented above are sufficient to compute the self-energy ([2| of small molecules [33]. As a test, we will 
compute the interacting Green's function by solving Dyson's equation ([lj. From this Green's function we can obtain 
the DOS and estimate the positions of the HOMO and LUMO levels. Here we illustrate this procedure in the case 
of benzene. This molecule has been chosen because extensive theoretical results and experimental data are available 
for it. Our calculations show a considerable improvement using the GW approximation as compared to the results 
obtained with plain DFT calculations using local or semi-local functionals. In general, for small molecules we find a 
reasonable agreement with experimental data and previous GW calculations of the ionization potentials and electron 
affinities. 

The input for our GW method has been obtained from calculations using the local density approximation (LDA) 
and the SIESTA package (TU] • SIESTA uses a basis set of strictly confined numerical atomic orbitals. The extension of 
these orbitals is consistently determined by an energy shift parameter. In general, the smaller the energy shift the larger 
the extension of the orbitals, although the procedure results in different cutoff radii for each multiplet of orbitals |35| . 
In the present calculations we have used the default double-^ polarized (DZP) basis, along with the Perdew-Zunger 
LDA exchange-correlation functional ZG\ and pseudo-potentials of the Troullier-Martins type [57]. Our calculations 
indicate (see table [I]) that it is necessary to use rather extended orbitals to obtain converged results for the HOMO 
and LUMO levels. For the most extended basis used here (determined from an energy shift of 3 meV) all the orbitals 
in benzene have a non-zero overlap and, in principle, the number of products of orbitals is 108(108+l)/2=5886. This 
number is reduced using the algorithm described in Ref. |12j . and the dominant product basis (see eq ( [16) ) only 
contains 2325 functions. The spectral functions have been discretized using a grid with — 1024 points in the 
range from —80 eV to 80 eV. The broadening constant has been set automatically to e — 1.5Ao> = 0.234375 eV. 
The frequency range was chosen manually by inspecting the non interacting absorption spectrum. The results of the 
calculation depend only weakly on the frequency range. 

Figure [2] shows the DOS calculated with different Green's functions. As one can see, the input Green's function 
Gq from a DFT-LDA calculation has a very small HOMO-LUMO gap. The Green's function G obtained with the 



instantaneous part of the self-energy (see eq 24) opens the HOMO-LUMO gap. This part of the self-energy S x 
incorporates the effect of exchange and is very important for small molecules. However, the gap is over-estimated 
as one can already anticipate from typical mean-field Hartree-Fock calculations. Correlation effects arc partially 
taken into account by the dynamical part of the GW self-energy. This brings the HOMO-LUMO gap closer to the 
experimental value. Our results stay also in agreement with other works using similar approximations (Go Wo on top 
of DFT-LDA) PUSH]. 

Apart from the GW approximations to the self-energy ([2]), our numerical method is controlled by precision param- 
eters of a more technical nature. Table [I] present the results for the ionization (IP) and electron affinity (EA) as a 
function of the extension of the atomic orbitals in the original LDA calculation as determined from the energy shift 
parameter [35 . An energy shift of 150 meV is usually sufficient to have an appropriate description of the ground-state 
properties of the molecules [TU] . However, we can see that our GW calculation requires more extended (smaller energy 
shift) orbitals. The slow convergence of the ionization potential of benzene with respect to the quality/completeness 
of the basis set has also been observed in the plane- wave calculations (using Wannier functions) of Ref. [Tj5] . 

Table [TJ also shows the results of calculations using one and two energy windows. The former calculation is more 
straightforward but requires the same density of frequency points at higher energies as in the region of interest near 
the HOMO and LUMO levels. The latter calculation uses two separate frequency grids: as described in section [V B| a 
lower resolution and a larger imaginary part of the energy are used for the whole spectral range, while high resolution 
and a small width are used in the low energy range to resolve the HOMO and LUMO levels. Thus, the second 
window technique requires the computation of both the response function and screened Coulomb interaction at far 
fewer frequencies than the one-window calculation. For instance, the one-window results presented above have been 
obtained with N u = 1024 frequencies, while the two-window results used only = 192 frequencies in both windows, 
implying a gain of a factor 2.7 in speed and in memory. The first and second window extend to 12.58 eV and 80 
eV, respectively. The first window is chosen as 2.5(£"dft lumo — -^dft homo)- The computational result depends 
only weakly on the extension of the first window. The second window is chosen manually as in the one-window 
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Binding energy, eV 

a) b) 

FIG. 2. a) Density of states of benzene computed from different Green's functions using as an input the results of a DFT-LDA 
calculation performed with the SIESTA package. A DZP basis set, with orbital radii determined using a value of the energy 
shift parameter of 3 meV, has been used. The results shown in this figure are obtained with a single energy window. GW X 
refers to the results obtained with only the instantaneous part of the self-energy (only exchange), while GW X c labels the results 
obtained with the whole self-energy (incorporating additional correlation effects), b) Ball and stick model of benzene produced 
with the XCrysDen package |39j . 



Energy-shift, 


One window 


Two windows 


meV 


IP, eV 


EA, eV 


IP, eV 


EA, eV 


150 


8.48 


-1.89 


8.48 


-2.01 


30 


8.71 


-1.45 


8.72 


-1.57 


3 


8.76 


-1.29 


8.78 


-1.41 


Experiment 


9.25 


-1.12 


9.25 


-1.12 



TABLE I. Ionization potentials and electron affinities for benzene versus the extension of the basis functions. The extension 
of the atomic orbitals is determined using the energy shift parameter of the SIESTA method [35]. Note that rather extended 
orbitals are necessary to achieve converged results. Differences associated with the use of the second window technique intro- 
duced in section ( |V B[ ) are of the order of 0.1 eV. The experimental ionization potential is taken from the NIST server |40j . 
The electron affinity of benzene is taken from [41] . 



calculation above. The broadening constant e has been set separately for each spectral window, using the default 
value e = 1.5Aw. The lower number of frequencies obviously accelerates the calculation and saves memory, while 
introducing very small inaccuracies in the low frequency region. According to table [I] the positions of HOMO and 
LUMO agree within 0.1 eV. Figure [3] shows that the second window leads to changes that are small, both in the 
HOMO and LUMO positions, and in the density of states at low energies. 

The calculations presented in this section needed a fairly large amount of random access memory (RAM). The 
amount of RAM increases as A? -2 with N the number of dominant products, which prohibits the treatment of larger 
molecules using the methods described above in a straightforward manner. However, as we will see in the next section, 
we can use a compression method that dramatically reduces the required memory. 



VII. COMPRESSION OF THE COULOMB INTERACTION 

As shown in Ref. [H] it is possible to solve the Petersilka-Gossmann- Gross equations |42 for time-dependent density 
functional theory (TDDFT) using a Lanczos type approach if, for example, we are only interested in the polarizability 
tensor of the system. In this way, we avoid keeping the entire linear response matrix x^iX^) i n the computer memory. 
Unfortunately, we were unable to find an analogous Lanczos type procedure for the self-energy matrix. However, 
we have found an alternative solution to this problem. It consists of taking into account the electron dynamics and 
keeping preferentially those dominant products that are necessary to describe x°„ in the relevant range of frequencies. 
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One window 
Two windows 




10 -5 

Binding energy, eV 



FIG. 3. The density of states of benzene computed with a uniformly discretized spectral function and using the second 
window technique. The peak positions are very weakly perturbed by using the two windows technique. The parameters of the 
calculation are identical with those of figure [2] The two windows technique allowed to reduce the number of frequency points 
from N u = 1024 to N u = 192. 



A. Defining a subspace within the space of products 



Consider the following closed form expression of the non interacting response °f ec l (191 



EF 



np — tie 



E.F 



ie-(E-F) 



V FF , where V FF = X*V*X?. 



(35) 



This is a well known expression, but rewritten in the basis of dominant products [13] . It must be emphasized that 
we do not use this equation to compute x^ui^) (it would require 0(N 4 ) operations) but this explicit representation 
of the exact non interacting response is nonetheless crucial for motivating our method of compression. 

Clearly x%(^) is built up from 0(N 2 ) vectors V? F . On the other hand, the entire space of products is, by 
construction, of only 0(N) dimensions. Therefore, there must be a significant amount of collinearity in the set of 
vectors V FF and a much smaller subset of such vectors should span the space where X°^( w ) acts. As candidates for 
the generators of this subspace, we sort the vectors {V FF } according to \E — F\ up to a certain rank iV ran k 



{A™} = subset of {V } limited according to \E — F\ < ^threshold) n = 1 . . . N, 



rank ■ 



(36) 



Here we treat {E, F} as electron-hole pairs, i.e. E < and F > 0. 

As a first test of whether the subspace carries enough information, we define a projector onto it 



ran vtn^uvvn. 
9 - X u v X v > 

P nl , = A™g m „A", where g r , 



(37) 



It can be shown without difficulty that Pf? is indeed a projector in the sense of P 2 
screened Coulomb interaction onto the subspace generated by the set {A"™} n= i.jv r8 



P. We can use it to project the 



(38) 



We must choose N tan k large enough so that the trace of the projected spectral density W / p r oj 0C tod( aJ ) ^ s sufficiently 
close to the original one. We checked that this works even for A ran k considerably smaller than the original dimension 
of the space of products. We can go further and reduce the dim ension of the subspace by eliminating collinear vectors 
from it. We do this by diagonalizing the matrix g mn in eq (37) and by defining new vectors 



(39) 
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To define the vector space {Z*}, we first discard the eigenvectors that correspond to eigenvalues A smaller than 
a threshold with respect to the Coulomb metric v^ u and we then normalize the remaining vectors, for simplicity. As 
a result of this procedure, we obtain a smaller set of vectors that we denote again by {Z nfl } with n = 1 . . . -/V su b ra nk> 
with iV su brank < N iax ± and the additional property that they are orthonormal with respect to the Coulomb metric 

Z^Z nfl = 5 mn , where Z£ = v^ v Z mv , for m,n = 1 . . . iV su brank- (40) 



Construction of the screened interaction from the action of the response function in the subspace 



From the preceding discussion we know that can be adequately represented in the previously constructed 
subspace {Z™} in the sense of 

Xfiv = Z m fiXmnZni>, with Xmn ~ ^mXfj,i/Z n - (41) 

To see which form the screened Coulomb interaction (|3| takes for such a density response Xui/> we write it as a series 



(42) 



Because x° a cts — by hypothesis — only in the subspace, the series may be simplified. Lets insert the representation 



of the response function of eq (41 ) into the series (42 1 



= xT + v^' [Z m ^x° mn Z nv ,] v»' v + v^' [Z mll ,x mn Z nv ,] v v '»" [Z m ^ X mn Zn„»] v v " v + 



then use the orthogonality property of the basis vectors Z mfl ( 40 ) and find 



W fiu nu , Z ti o r c , o , o o , 

vy ' ZJ mA.mr |_ rn ^ Am ~ ArsAsn ' 



(43) 



RPA r?V 

mAmn n' 



Z^Y 
m A 



Here we introduced the new response function Xmn A = {&mk — Xmk) Xkn- From the preceding arguments we 
conclude that the dynamically screened Coulomb interaction W^ v can be computed in terms of the response function 
Xm« A within the previously constructed subspace and the matrix inversion in this smaller space is of course much 
cheaper than in the original space. This is a welcome feature — the number of operations for matrix inversion scales 
with the cube of the dimension and a compression by a factor 10 will lead to a 1000 fold acceleration of this part of 
the computation. 

It is important to note that, although an energy cut-off ^threshold is use d to choose the relevant {V® F } vectors, 
high frequencies components of the response and the screened interaction are explicitly calculated, -^threshold only 
serves to construct the frequency independent basis vectors {Z£} according to eq (39 1. This basis is later used to 
calculate the response function x°„(u;) in the whole frequency range (see eq (41 )). Of course, we can expect that, 



if -^threshold is chosen too small, the ability of the compressed basis to represent the high energy components of the 
response will eventually deteriorate. However, we are interested in the low energy excitations of the system and, as 
we will show in the next subsection, those can be accurately described using values of -^threshold that allow for a large 
reduction in the size of the product basis. Furthermore, it is also important to note that the instantaneous self-energy 
S x , for which a compression criterion based on our definition of -EWcshoid is dubious, is calculated within the original 
dominant product basis, i.e. before this non local compression is performed. 



C. The compression in the case of benzene 



The non local compression depends on two parameters: i) the maximum energy -^threshold of the Kohn-Sham 
electron-hole pairs {V® } in eq (36) , and ii) the eigenvalue threshold A for identifying the important basis vectors 
{Zi;\ in eq Q. 

Table |H| shows the electron affinity of benzene as a function of ^threshold & n d A. The computational parameters have 
been chosen as in section [Vl] and the energy shift to define the extension of the orbitals is 3 meV. Table |H| illustrates a 
general feature that we have found in many test for several systems: N TaB k can be chosen of the order of the number 
of atomic orbitals iV or b. We have found that iV ran k ~ 5iV or b usually guarantees a converged result for the HOMO and 
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A = 1Q-' 2 


A = 10" a 


A = 10~ 4 


^threshold — 10 eV 
^threshold = 20 eV 
^threshold = 40 eV 


2.48 (33) 
1.38 (96) 
1.41 (132) 


2.46 (37) 
1.39 (133) 
1.41 (192) 


2.45 (39) 

1.40 (171) 

1.41 (279) 



TABLE II. Electron affinities for benzene versus the compression parameters ^threshold and eigenvalues cutoff A in eq (39 1. 
In brackets the dimension of compressed subspace is given JV su b ran k- The dimension N la ,nk is governed by ^threshold an d was 
39, 297 and 765 for -^threshold = 10 eV, i?thr CS hoid = 20 eV and -^threshold = 40 eV, accordingly. Number of atomic orbitals 
A'orb = 108, while the number dominant products is jVp ro d = 2325. 



LUMO levels. In any case, the number of relevant linear combinations iV su b ra nk was always much smaller than the 
number of dominant functions, with a typical compression ratio of ten or more. 

To further illustrate the quality of the basis, we will compare the trace of the original screened interaction with 
the trace of the projected screened interaction for the benzene molecule. The result of the comparison can be seen 
in figure [4] In this test calculation, the dominant product basis consists of 921 functions, while the compressed basis 
contains only 248 functions. Examples of compression for larger molecules will be presented in section [Xj 
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FIG. 4. a) Comparison of the screened interaction calculated for be nzen e using our original dominant product basis and the 
screened interaction projected to a compressed product basis (see eq (38 1). We plot the sum of all the matrix elements of the 
imaginary part of the screened interaction, b) A plot of the difference of the functions represented in panel a. The change 
in spectral weight of the screened interaction due to compressing the space of dominant products is seen to be small. Please 
notice the different scales of the y-axis in both panels. 



The examples presented in this section show that the screened Coulomb interaction can be effectively compressed. A 
practical algorithm that uses the non local compression and maintains the 0(N 3 ) complexity scaling of the calculation 
will be presented in the next section. 



VIII. MAINTAINING 0(/V 3 ) COMPLEXITY SCALING BY COMPRESSING / DECOMPRESSING 



The favorable 0{N 2 ) scaling of the construction of the uncompressed non interacting response x°^i> is due to its 
locality. On the other hand, we need compression for ^ to fit into the computer memory and the compressed Xmn * s 
no longer local. To satisfy the two mutually antagonistic criteria of i) locality (for computational speed) and ii) small 
dimension (to fit into the computer memory) we shuttle back and forth as needed between the uncompressed/local 
and the compressed/non local representations of the response x° and of the screened interaction W. Both compression 
and decompression are matrix operations that scale as 0(N 3 ) and this, along with the matrix inversion in eq (43) 
in the computation of the screened Coulomb interaction, and the computation of the spectral densities ptu(s) is the 
reason why our implementation of GW scales as 0(N 3 ). 
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A. A construction of the subspace response in 0(N 3 ) operations 
Let us describe an efficient construction of the response Xav and its compressed counterpart Xmn that, besides, 



gives us an opportunity to describe our use of frequency and time domains during the calculation. Consider eq (32) 
that involves convolutions of the spectral functions pt c (u>) and p~ d {u) = p~ d (—uj). To make use of the convolution 
theorem, we will first compute the spectral function of the non interacting response a fll ,(s) in the time domain 



2tt ' J a " r " cv 2tt ./n v aax ' 2tt 

ab ^+ /j.\jrcd- 



= ^pt(t)V^p- d (t). (44) 

In other words, we prepare the use of the FFT driven convolution by first computing the Fourier transforms of the 
electronic spectral densities p ± and once a^(t) is computed, we return to a llv (s) by an inverse Fourier transformation. 
This is nothing else but the fast convolution method with the Fourier transform of the spectral densities p^(w) and 



p ad {u) carried out prior to the tensor operations in eq (44) 



Above we saw how to compute the spectral function of the non interacting density response. However, as we 
mentioned before, we cannot easily store this information in the memory of the computer and we must therefore 
compress this quantity as soon as it is found to avoid over flooding the computer memory. An efficient way to do this 
is to compute a^vit), the spectral function of the non interacting response in the time domain, in a "time by time" 
fashion, with the time variable t in the outer loop. Fortunately, the compression of the response Xuv according to 



equation (41) can be done on the level of its spectral function a MI ,(i) separately for each time t. 



B. A construction of the self-energy in 0(N 3 ) operations 



Although we use the spectral function given by eq (30) to compute the self-energy in the second of eqs (29), it 
is useful to think also of eq (21) that corresponds to the Feynman diagram of Figure [l] and which has the same 
locality properties. Please recall that the product vertex in eq (16) is sparse and local and that the indices {a,a',p} 
and {6,6',^} must each reside on a single pair of overlapping atoms. Once the indices a, b of the self-energy are 
specified, there are only O(N ) possibilities of choosing the remaining indices. Therefore, the calculation of £ ab (£) 
requires asymptotically 0(N 2 ) operations provided that the screened Coulomb interaction W^ u in a basis of localized 
functions is known. However, the local screened Coulomb interaction in the original space of dominant products 
does not fit into the computer memory as opposed to the compressed, but non local, response Xm« A that we store 
(see eq (43)). We may, however, regain locality by decompressing Xmn A at the cost of 0(N 3 ) operations, using the 



identity W^y naroical = Z m x m F n AZ n in eq Q. As we cannot keep W^ namical in the computer memory, we must try to 



"decompress Xmn on the nv "- To do this, let us transform the first of eqs (30) into the time domain. For instance, 
for the positive part of the spectral density cr+ b (i) of the self-energy, we find 



*?(t) = 2nV™'p+ b ,(t)V v b ' b 1 i?(t). 

Again, the representation in time of p^, b , (t) is prepared only once. However, the transform j^_ u (t) — — \ Z^Imx^^ (t) Z\ 
for all times does not fit into the computer memory. Therefore we also decompress "/^"(t) time by time by letting the 
time t run in the outer loop, by computing 7+"(t) via decompression for a single time, and by storing only the result 
c±{t) for each time. Once we have computed cr±{t) for all times, we can find <J±(s) from it. 



IX. A SUMMARY OF THE COMPLETE ALGORITHM 



At this point, it is useful to briefly recapitulate the different steps of our implementation of Hedin's GW approxi- 
mation. It consists of the following steps: 

1. Export the results of a DFT code that uses numerical local atomic orbitals as a basis set. Here we use the 
SIESTA code |T0 ], but other codes like the FHI-AIMS code [44] could also be used. 

2. Set up a basis of dominant products in O(N) operations. Here we use the method of Ref. [1"2"] . 
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Naphthalene 


Anthracene 


Energy-shift, 


One window 


Two windows 


Two windows 


meV 


IP, eV 


EA, eV 


IP, eV 


EA, eV 


IP, eV 


EA, eV 


200 


7.24 


-0.68 


7.27 


-0.79 


6.44 


0.20 


20 


7.61 


-0.083 


7.67 


-0.18 


6.89 


0.77 


Experiment 


8.14 


-0.191 


8.14 


-0.191 


7.439 


0.530 



TABLE III. Ionization potentials and electron affinities for naphthalene and anthracene and their dependence on the extension 
of the atomic orbitals. For naphthalene we compared the results obtained with spectral functions discretized in one or two 
windows. The experimental data has been taken from the NIST server [40]. For naphthalene and anthracene, vertical ionization 
potentials are not available at the NIST database. Therefore we give experimental ionization energies including effects of 
geometry relaxation. 



3. Set up a space of reduced dimension where the screened Coulomb interaction will act and exploiting the low 
effective rank of this set. Such a subspace is determined by a set of iV ran k vectors V^ F that correspond to 
electron- hole pairs with a predetermined maximum value of \E — F\. Further compression is obtained by 
diagonalizing the Coulomb metric projected onto this subspace. This step requires 0(N 3 ) operations. 

4. Choose low and high energy spectral windows and a frequency grid. Prepare the electronic spectral density 
Pab(s) in these two windows from the output of the DFT calculation. 

5. Find by constructing and compressing the local "on the fly" in 0(N 3 ) operations and by solving 
for Xmn A f° r au frequencies in 0(N 3 ) operations. The construction must be done in two frequency windows. 
Truncate the spectral data where needed in order to avoid double counting and store in the two windows. 

6. Find the spectral function of the self-energy by decompressing Xmn A ~^ W^" "on the fly" in 0(N 3 ) operations 
and by convolving it with the electronic spectral function. Again this must be done in two frequency windows 
and the results must be combined consistently. 

7. Construct the self-energy from its spectral representation. 

8. Solve Dyson's equation and find the density of states from the interacting Green's function. Obtain the desired 
spectroscopic information from the density of states. 

Results obtained with the above algorithm will be discussed in the next section. 



X. TESTS FOR MOLECULES OF INTERMEDIATE SIZE 



The compression technique has been carefully tested in the case of the benzene molecule. The tests show excellent 
agreement of the density of states computed with and without compression. In this section, we will consider larger 
molecules such as the hydrocarbons naphthalene and anthracene |34) . These molecules are well known to differ in 
their character as electron acceptors: naphthalene, like benzene, has a negative electron affinity, while anthracene is 
an electron acceptor with positive electron affinity. 

A compression of the dominant product basis is necessary to treat the molecules considered in this section. These 
molecules are too large for a calculation without compression on ordinary desktop machines because of memory 
requirements. For naphthalene the dominant product basis contained 4003 functions, which were reduced to 433 
functions after compression. In the case of anthracene, the dominant product basis contained 5796 functions, while 
the compressed basis had only 598 functions. 

Table III shows our results for naphthalene and anthracene. The computational details were similar to those used 
for benzene and already described in section |VI| The two- window results were obtained with frequency grids of only 
128 points for naphthalene (in the ranges ± 8.32 eV and ± 80 eV) and yet it provides an accuracy on the 0.1 eV 
level, while the one- window calculation is done again with 1024 frequencies (in the range of ± 80 eV). In the case of 
anthracene, results using frequency grids of 256 points (in the ranges ± 16 eV and ± 90 eV) are presented. 

Again we find a large improvement over the position of the Kohn-Sham levels in a DFT-LDA calculation. The 
agreement with the experimental data is certainly improved, although there are still significant deviations, particularly 
with respect to the reported ionization potentials. Interestingly, however, our calculations recover the important 
qualitative feature of anthracene being an electron acceptor. In the case of anthracene, the results obtained with 
our most extended basis orbitals (energy shift of 20 meV) are in excellent agreement with the recent calculations by 
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FIG. 5. a) Density of states for naphthalene. The results have been obtained with our most extended basis orbitals (corre- 
sponding to an energy shift of 20 meV [35]). We can appreciate the accuracy of the second window technique, b) Ball and 
stick model of naphthalene produced with the XCrysDen package [39] . 



Blase et al. [5j . In the case of naphthalene we can see that the two frequency windows technique introduces only tiny 
differences (below 0.1 eV) in the positions of the HOMO and LUMO levels. 

The corresponding DOS is shown in figures [5] and [6j One can see in figure [6] that it is the dynamical part of the 
self-energy, including correlation effects, that turns our theoretical anthracene into an acceptor, while including only 
the instantaneous self-energy predicts anthracene to be a donor. 




Binding energy, eV 

a) b) 

FIG. 6. a) Density of states for anthracene. The results have been obtained using the extended basis orbitals corresponding 
to an energy shift of 20 meV |35| . Here we compare calculations using the instantaneous (exchange-only) self-energy and 
the full self-energy (including correlation effects). The correlation component of the self-energy is crucial to reproduce the 
experimental observation that anthracene is an acceptor. In contrast, the exchange-only calculation locates the LUMO level 
above the vacuum level, b) Ball and stick model of anthracene produced with the XCrysDen package |39| . 

These results for molecules of modest size are just a first application of our algorithm. With its favorable scaling, 
our method aims at GW calculations for larger molecules of the type used in organic semiconductors. However, before 
carrying out such studies, we should reduce the initial number of dominant products i. e. before any compression is 
applied to it. 



XI. CONCLUSIONS AND OUTLOOK 

In the present paper we have described our approach to Hedin's GW approximation for finite systems. This 
approach provides results for densities of states and gaps that are in reasonable agreement with experiment and it 
requires only modest computer resources [34 for the systems presented here. The complexity of our algorithm scales 
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asymptotically as the third power of the number of atoms, while the needed memory grows with the second power of 
the number of atoms. We hope that these features, along with a further reduction of the size of the basis describing 
the products of localized orbitals, will allow to apply our method to describe the electronic structure of large molecules 
and contribute to an ab-initio design of organic semiconductors for technological applications. 

The algorithm described here is built upon the LCAO technique |10j and uses a previously constructed basis 
in the space of orbital products that preserves their locality and avoids fitting procedures |12j . Moreover, a (non 
local) compression technique has been developed to reduce the size of this basis. This allows to store the whole 
matrix representation of the screened Coulomb interaction at all time/frequencies in random access memory while 
significantly reducing the computational time. The time (and frequency) dependence of observables is treated with 
the help of spectral functions. This avoids analytical continuations and allows for operations to be accelerated by 
the use of FFTs. As a useful byproduct of our focus on spectral functions we obtain, as primary result, an electronic 
spectral function of the type observed in photo-emission and from which we then read off the HOMO and LUMO 
levels. 

We have applied our method to benzene, naphthalene and anthracene. As expected, we find that our estimations 
of the HOMO and LUMO positions and the corresponding gaps are significantly improved over the results obtained 
from the Kohn-Sham eigenvalues in a plain DFT-LDA calculation. Our results approach the experimental data but, 
as observed by other authors these "single-shot" GW-LDA calculations (or GoWo-LDA using a more standard 
terminology) still present sizeable deviations from the measured ionization potentials and electron affinities. In general, 
our results are in good agreement with previous GoWo-LDA calculations for similar systems [3[T9l|38]. Thus, we 
expect further improvements by iterating our procedure until self consistency or, as suggested by other authors in the 
case relatively small molecules [5] [201 EH] j by using Hartree-Fock results as an input for our Go Wo calculations. For 
periodic systems it is well known that GoWo-LDA systematically underestimates the size of the gaps of semiconductors. 
The best results so far were found using the so-called "improved quasi-particle method" [21 [46] . A realization of this 
method in our framework should also improve the precision of our results. 

The method presented in this paper depends crucially on the quality and size of the original LCAO basis. A possible 
limitation is that the typical LCAO basis used in electronic structure calculations are constructed and optimized in 
order to describe ground-state properties [35] . However, it is possible to optimize an LCAO basis, for example using 
a technique similar to that described in Ref. [17], to represent electronically excited states. This will increase the 
accuracy and applicability of the method and could even allow to reduce the size of the original LCAO basis used to 
represent the electronic states. Moreover, by comparing our basis with that of other authors, there are indications 
that the (local) basis of dominant products used in this paper can be reduced in dimension without changing the 
physical results [5]. Such a reduction should lead to an important improvement of the prefactor in our implementation 
of GW, but also, as a side effect, introduce a similar acceleration in our published TDDFT algorithm [T3] that is 
already competitive, in its present form, with other TDDFT codes. 

The quantities calculated in the presented algorithm can be useful in other branches of many-body perturbation 
theory. For instance, the screened Coulomb interaction is a crucial ingredient of the Bethe-Salpeter technique that is 
needed to study excitons and the optical response of excitonic systems. In this context it is interesting to note [15] 
that the solution of the Bethe-Salpeter equation scales as 0(N 3 ) for clusters of size N, at least when suppressing the 
dynamic part of the fermion self-energy and the dynamic part of the screening of the Coulomb interaction. Calculations 
of the transport properties of molecular junctions |49j are another possible application of the GW approach described 
here. 
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